Aggregate stressors to 10 km Mollweide
Compiled on Tue Jun 16 14:37:47 2020 by ohara
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
source(here('common_fxns.R'))
dir_str_mol <- '/home/shares/ohi/git-annex/impact_acceleration/stressors'
dir_str_10km <- file.path(dir_bd_anx, 'layers/stressors_10km')1 Summary
Collect stressor layers from CHI 2019 and reproject/aggregate them to the native resolution of the species range maps. This is a simple aggregation 11x as is done for the ocean, EEZ, MEOW, etc rasters.
1.1 Loop over all stressor layers, aggregate to spp CRS
Read in each raster, aggregate to a raster at the species projection and resolution. Because the projection is Mollweide for the source and target, and the base raster for species maps was created the same way, we can simply aggregate the original Mollweide CHI maps up by the same factor. For spp maps, the aggregation factor was 11\(\times\) to approximate 100 km2 cells: \(.934 \times 11 = 10.28 \text{ km}; \; (.934 \times 11)^2 = 105.6633 \text{ km}^2\).
Notes:
- here we are aggregating using mean value.
- land-based stressors will be masked to the coastal 3nmi buffer before aggregating. Real zeros in coastal zones will still be counted, but zeros farther offshore will be ignored.
- for sea surface temperature, we will mask out sea ice prior to aggregating. The original file for the sea ice mask comes from:
mazu:git-annex/Global/NCEAS-Pressures-Summaries_frazier2013/ice_mask_resampled(both.griand.grdfiles needed)
### List directories that end in "final"
str_dirs <- list.files(dir_str_mol, recursive = TRUE, full.names = TRUE) %>%
dirname() %>%
unique() %>%
.[str_detect(., 'final') & !str_detect(., 'archive')]
# [1] "archive/art_fish/final" "archive/comm_fish/final/dem_dest"
# [3] "archive/comm_fish/final/dem_nondest_hb" "archive/comm_fish/final/dem_nondest_lb"
# [5] "archive/comm_fish/final/pel_hb" "archive/comm_fish/final/pel_lb"
# [7] "archive/sst_old/final" "art_fish/final"
# [9] "benthic_structures/final" "comm_fish/final/dem_dest"
# [11] "comm_fish/final/dem_nondest_hb" "comm_fish/final/dem_nondest_lb"
# [13] "comm_fish/final/pel_hb" "comm_fish/final/pel_lb"
# [15] "direct_human/final" "land_based/final/nutrient"
# [17] "land_based/final/organic" "light_pollution/final"
# [19] "oa/final" "shipping/final"
# [21] "slr/final" "sst/final"
# [23] "uv/final" ### identify land-based stressors, and load coastal mask and sea-ice raster.
land_based_str <- get_str_cats() %>%
filter(category == 'land-based') %>%
.$stressor
coastal_str <- c('slr', land_based_str)
coastal_mask <- raster(file.path(dir_bd_anx, 'spatial/three_nm_offshore_mol.tif'))
seaice_mask <- raster(file.path(dir_bd_anx, 'spatial/sea_ice_mask_1km.tif'))
for(dir in str_dirs) {
### dir <- str_dirs[14] ### slr
### dir <- str_dirs[13] ### shipping
message('Processing stressors in ', dir)
files <- list.files(dir, pattern = '.tif$', full.names = TRUE)
# for(f in files) {
tmp <- parallel::mclapply(files, mc.cores = 12,
FUN = function(f) {
### f <- files[1]
str <- str_replace_all(basename(f), '_[0-9].+', '')
str_10km_file <- file.path(dir_str_10km,
str_replace(basename(f), 'mol', 'mol10km'))
if(!file.exists(str_10km_file)) {
message('Processing ', basename(f))
str_rast_orig <- raster(f)
### check for coastal or sst stressors
if(str %in% coastal_str) {
message('Masking ', basename(f), ' to coastal only')
str_rast_orig <- mask(str_rast_orig, coastal_mask)
}
if(str == 'sst') {
message('Masking ', basename(f), ' to remove sea ice')
str_rast_orig <- mask(str_rast_orig, seaice_mask)
}
str_rast_10km <- raster::aggregate(str_rast_orig,
fact = 11, fun = mean,
filename = str_10km_file,
overwrite = TRUE)
### check to make sure results match spp rasts
# ocean_area_rast <- raster(here('_spatial/ocean_area_mol.tif'))
# raster::compareRaster(str_rast_10km, ocean_area_rast,
# extent = TRUE, rowcol = TRUE, crs = TRUE,
# res = TRUE, orig = TRUE)
### TRUE!
} else {
# message('File exists: ', str_10km_file, '... skipping!')
}
})
}1.2 Map all stressors in gifs
library(animation)
make_gifs <- function(map_stack, filename, layer_names = NULL) {
if(is.null(layer_names)) layer_names = names(map_stack)
x <- values(map_stack)[values(map_stack) > 0] %>%
min(na.rm = TRUE)
power_adj <- -1 * round(log10(mean(values(map_stack), na.rm = TRUE) * 2)) + 1
message('Adjusting by power of ', 1/power_adj)
map_stack_adj <- map_stack^(1 / power_adj)
ramp_breaks <- c(0, seq(x, 1, length.out = 100))
ramp_labels <- seq(0, 1, .2)
ramp_colors <- c('grey40', hcl.colors(n = 100, palette = 'viridis'))
capture.output({
saveGIF({
for(i in 1:nlayers(map_stack_adj)){
plot(map_stack_adj[[i]],
col = ramp_colors,
breaks = ramp_breaks,
axes = FALSE,
axis.args = list(at = ramp_labels^(1 / power_adj), labels = ramp_labels),
main = layer_names[i])
}},
interval = 0.5, movie.name = filename,
ani.width = 700, ani.height = 420)
})
return(invisible(NULL))
}str_10km_files <- list.files(file.path(dir_bd_anx, 'layers/stressors_10km'),
full.names = TRUE)
years <- 2003:2013
str_map_df <- data.frame(str_file = str_10km_files) %>%
mutate(str = str_replace(basename(str_file), '_[0-9].+', ''),
year = str_extract(basename(str_file), '[0-9]{4}') %>%
as.integer()) %>%
group_by(str) %>%
filter(min(year) <= 2003 & max(year) >= 2013) %>%
ungroup() %>%
filter(year %in% years)
stressors <- str_map_df$str %>% unique()
for(stressor in stressors) { ### stressor <- stressors[14]
### Animate the results
gif_file <- here('_setup', sprintf('figs/str_movie_%s.gif', stressor))
str_maps <- str_map_df %>%
filter(str == stressor) %>%
.$str_file
if(!file.exists(gif_file)) {
message('creating ', gif_file)
rast_files <- str_maps
map_stack <- stack(rast_files)
make_gifs(map_stack,
filename = gif_file,
layer_names = paste(stressor, years))
}
cat(sprintf('', gif_file))
}
